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For a distinguished water molecule, the solute water, we assess the contribution of each coordina- 
tion state to its excess chemical potential, /x™, using a molecular aufbau approach. In this approach, 
we define a coordination sphere, the inner-shell, and separate /x™ into packing, outer-shell, and local 
chemical contributions; the coordination state is defined by the number of solvent water molecules 
within the coordination sphere. The packing term accounts for the free energy of creating a solute- 
free coordination sphere in the liquid. The outer-shell term accounts for the interaction of the solute 
with the fluid outside the coordination sphere and it is accurately described by a Gaussian model 
of hydration for coordination radii greater than the minimum of the oxygen-oxygen pair correlation 
function. Consistent with the conventional radial cut-off used for defining hydrogen-bonds in liq- 
uid water, theory helps identify a chemically meaningful coordination radius. The local chemical 
contribution is recast as a sum over coordination states. The n"^ term in this sum is given by 
the probability of observing n water molecules inside the coordination sphere in the absence of the 
solute water times a factor accounting for the interaction of the solute with the inner-shell solvent 
water molecules. Using this molecular aufbau expansion, we monitor the change in the chemical 
contribution due to the incremental increase in n. We find that though four water molecules are 
needed to fully account for the chemical term, the first added water, the n = 1 coordination state, 
accounts for nearly half the chemical term. Our results emphasize the need to acknowledge the 
intrinsic occupancy of a solute-free coordination sphere together with solute-solvent interactions in 
rationalizing the tetrahedral coordination of the solute water. 



I. INTRODUCTION 

Computer simulation studies of water have signifi- 
cantly contributed to our understanding of the struc- 
ture of the liquicP^. In broad agreement with X-ray 
and neutron scattering experiments, these studies show 
that at standard conditions of temperature and pressure 
each water molecule in the liquid donates approximately 
two hydrogen bonds. But in contrast to ice, the hydrogen 
bonding is characterized by substantial disorder^. For ex- 
ample, as Fig.[l]shows for two common simulation models 
of water, the most probable coordination state of a wa- 
ter molecule is four, but with some variability about this 
coordination state. 

Recent studiePI^ based on X-ray absorption and emis- 
sion spectroscopy, however, suggest that the majority of 
water molecules in the liquid donate only one, strong hy- 
drogen bond. Since the well-studied simulation models 
of water (Fig. [ij describe reasonably welP^l^l the phase 
behavior of the fluid at normal temperature and pres- 
sure, this suggestion naturally raises the question of how 
the coordination states of water are related to the ex- 
cess chemical potential, a basic descriptor of the phase 
behavior. However, even for the well studied simulation 
models, the contribution of the coordination states to the 
free energy is not known. In this article, we present the 
theoretical framework to address this issue and apply it 
to the SPC/E'* water model. 

The theoretical framework anchoring the present study 
derives from a quasi-chemical organization of the poten- 
tial distribution theorenPiltS!^ Within this approach a 
coordination region is defined around the solute to sort 
local solvent configurations from bulk arrangements of 
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FIG. 1. Probability of observing various number n of wa- 
ter molecules within a sphere of radius 3.3 A around a distin- 
guished water molecule for two common simulation models of 
water, SPC/F''^and TIP4P/2005^. The position of the oxygen 
atoms defines the position of the water molecules. The mini- 
mum of the oxygen-oxygen pair-correlation function defining 
the first hydration shell is at 3.3 A. 

the liquid. Implicit in this sorting is the separation of 
the strong and specific local interactions from the non- 
specific solute-bulk interaction^^^"^^'. The excess chem- 
ical potential of the solute is thus given by the sum of: 
(1) a local packing contribution accounting for the work 
of creating a cavity in the liquid, (2) chemical contribu- 
tion accounting for the specific local interactions within 
the coordination region, and (3) contributions from the 
long-range interactions. 
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The local contribution can be expressed as a sum over 
coordination statej^^ ^-*-^!--^, where a coordination state 
is labelled by the number n of solvent molecules in the co- 
ordination region. This allows one to examine the effect 
of progressively adding solvent molecules to the solute's 
coordination region while monitoring its excess chemical 
potential. Using this molecular aufbau approaclP^, ear- 
lier we studied the hydration of some monovalent ions. 
For chemically reasonable coordination radii, the molec- 
ular aufbau approach showed that the local chemical 
contribution is almost entirely obtained by coordination 
states below the most probable coordination state. Fur- 
ther, only a subset of those states make a dominant con- 
tribution to the thermodynamics of hydration. For exam- 
ple, the most probable coordination number of K"'"(aq) 
is n = 6, but almost all of the local chemical contri- 
bution is obtained by the n — A coordination state. The 
most probable coordination state of a distinguished water 
molecule is n = 4 for SPC/E water (Fig.[l]). The earlier 
study then suggests that n < 4 states are important in 
accounting for specific interactions, but the contribution 
of each individual state remains to be determined. 

The rest of the article is organized as follows. In the 
following section we summarize the statistical thermo- 
dynamic theory. Since coordination states with small n, 
especially the state n = 0, may not be accessible via ther- 
mal fluctuations alone, we develop an ensemble reweigh- 
ing technique to uncover these low-coordination states. 
We study coordination radii up to 5.5 A. For a coordina- 
tion radius of 3.5 A, slightly greater than the first mini- 
mum of the oxygen-oxygen pair correlation function, the 
long range contribution is well-described by a Gaussian 
model of hydratiorf^SHIIl that is, the outer-shell contribu- 
tion is given by non-specific interactions. For this radius, 
the local chemical contribution is entirely accounted for 
by the n — A coordinate state, but the n — 1 state makes 
a contribution that is twice as large as the n = 2 state. 
Thus even for a fluid characterized as tetrahedral, we find 
that the four coordinating water molecules are different 
in their relevance to the thermodynamics. 

For a coordination radius of 5.5 A, none of the coordi- 
nation sphere water molecules directly bond with the dis- 
tinguished water molecule. They are instead pulled closer 
to the interface with the bulk. Investigating this aspect 
reveals a central insight of this paper: the occupancy 
characteristics of a solute-free coordination sphere are an 
important factor in the observed coordination structure 
of a solute-water. 



II. THEORY 

The devel opment below is anchored by the following 
identitji™ 

xn = p„e-'5[^" (")-'^w 1 , (1) 

where a:„ is probability of observing n solvent molecules 
within a defined coordination sphere, the inner shell. 



around the solute w. Here the solvent is water, and the 
position of the water molecule is defined by the position 
of the water oxygen atom. Further, in the present study, 
the solute w(= II2O) is also a water molecule. The solute 
water is termed the distinguished water molecule and the 
solvent molecules within the coordination sphere will be 
termed the solvent ligands. 

Pn is the probability of observing n solvent ligands 
within the inner shell, but with the solute-solvent in- 
teractions turned-off. /i^ is the solute's excess chemi- 
cal potential at temperature T (/3~^ = k^T). is 
the excess chemical potential conditioned on there be- 
ing only n solvent ligands within the inner-shell. Thus 
Pn, the intrinsic propensity of solvent ligands to occupy 
the inner-shell volume, is modified by solute-solvent in- 
teractions, defined by — /x^, to give the observed 
coordination number distribution Xn- 

The n = instance of Eq. [T] 

= kBTlnxo - kBThipn + fil^{0) (2) 

has a particularly clear physical meaning. The excess 
chemical potential of the solute, fj,'^, is given by the sum 
of three contributions. A packing term, — fceTlnpo? 
counting for the fr ee en ergy of creating an empty inner 
shell in the solvenlP2121l. (For convenience, we consider 
a spherical inner-shell of radius A.) Next the solute is 
placed in this cavity; the contribution due to its inter- 
action with the medium outside the inner shell is given 
by /z™(0). The absence of close solute-solvent contacts 
can make /i5f(0) particularly easy to model; for suffi- 
ciently large A, the interaction of the solute with the 
outer-shell solvent is entirely non-specific and a Gaussian 
model of the distribution of binding energies suffices ^'^ ^3. 
fceT In xo is the free energy gained by allowing solvent lig- 
ands to enter the inner shell. This contribution accounts 
for all the local, specific solute-solvent interactions. Note 
that each individual contribution on the right-hand side 
of Eq. [2] depends on A, but the left-hand side does not. 

Combining Eqs. [T] and [2] gives the molecular aufbau 
expansion of the inner-shell chemical contributiorP^ 

In = - In y ^e-'3[^w (")-A^^-(0)] . (3) 

Wn = ~ (0) is the free energy to compose a 

solute plus n-ligand cluster within the inner-shell and is 
a measure of local solute-solvent interaction. Each suc- 
cessive addition of a solvent ligand to the inner shell — a 
building up of solvent around the solute — leads to the 
inclusion of an additional contribution in the sum above. 
The contribution due to the addition of the i^^ solvent 
molecule {i > 1) is thus given by Inxo^i — Inxo.i-i, where 
Ina^o.i is the i**^ partial sum. 

Eqs. [I] and [3] make it clear that we need the entire {pn} 
and {x^ distributions to use the aufbau approach to as- 
sess the importance of individual coordination states. 
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A. Reweighing technique and occupancy statistics 



Neat Water 



For small coordination radii, {pn} and can be 

obtained directly from the simulation record or by us- 
ing a maximum entropy appr oach w ith known mean and 
variances of the distribution^^HHl or by stratified sam- 
pling approacheJ^. However, as the coordination radius 
increases, the lower coordination states are no longer ac- 
cessed and non-Gaussian behavior of {p„} and {xn} lim- 
its the use of the maximum entropy approaclP^. Inspired 
by the method of expanded enscmbles^^^^, we develop a 
reweighed sampling approach to construct the {pn} and 
{xn} distributions. (The method is essentially an impor- 
tance sampling approach.) 

Let Q{N,V,T,j) represent the canonical partition 
function when the solute has precisely j of the N total 
solvent molecules in the inner shell. Then, 



To determine {a;„}, {p„} for coordination spheres of 
different sizes and to calculate the van der Waals contri- 
bution to hydration, we simulated neat water for 6x10^ 
sweeps, where each sweep consists of an attempted move 
for every water molecule. The first 3x10^ sweeps were 
used for equilibration, during which the maximum al- 
lowed linear displacement and angular rotation were op- 
timized to yield an acceptance ratio of 0.3. For the next 
3x 10^ sweeps the optimized values of linear displacement 
and angular rotations were used unchanged. Configura- 
tions were saved every 10 sweeps for analysis. Except for 
the number of sweeps, a similar strategy was followed for 
the reweighed simulations as well. 



Pj = 



QiN,V,T,j) 
j:^QiN,V,T,j) 



(4) B. Reference values of hydration free energy and 

outer term 



Now if an energetic penalty rjj is applied to coordi- 
nation state j and the system is resampled, the revised 
probability of occurrence of state j is 



Pj 



Q{N,V,T,j) 



Q{N, V, T, j) E, QiN, V, T, j)ePv, ' 



(5) 



where the overbar denotes quantities in the reweighed 
sample. Clearly, for any two states j and r, we have 



Pi 

Pr 



Pj 



Pre 



(6) 



Choosing r to be the most probable state ensures that 
Pr is well-determined from the simulations. We can then 
adjust rjj relative to rjr to amplify the occurrence of the 
j state at the expense of the r state. A similar approach, 
but with a different set of weights {?7j\w}, can be used to 
obtain {a;„}, the probability distribution in the presence 
of the solute. 



III. METHODS 

We study liquid water under N VT c onditions using 
Metropolis Monte Carlo simulations'^SEIl. The cubic sim- 
ulation cell comprises 901 water molecules at the number 
density 33.33 nm"^. Water is modeled with the SPC/E 
potential. Simulations involving a distinguished water 
molecule were performed by fixing one of the 901 waters 
at the centre. Electrostatic interactions were modeled by 
the generalized reaction field (GRF^^^Hal] approach with 
a screening length equal to half the box-length. Lennard- 
Jones interactions were truncated at half the box-length. 
Statistical uncertainties in data are at the level of a k^T. 



To establish a reference for the hydration free energy 
of a water, we followed our earlier work^*. Electrostatic 
contributions to the excess chemical potential were ob- 
tained using a 2-point Gauss-Legendre quadratur^^, and 
the van der Waals contri bution was obtained using the 
histogram overlap methocP^'^SI Pqj- \ > 3 5 the van 
der Waals contribution to M^(0) was obtained using a 
Gaussian modeP^I^. A similar strategy was adopted 
in calculating the hydration free energies A*^(ra) for all 
states occurring within a 4.0 A coordination sphere. 



C. Estimation of {pn} and{a;„} 

For constructing {p„}, a cubic grid consisting of 343 
points was defined within the simulation cell. Observa- 
tion spheres of various radii were defined around each lat- 
tice point and the population of water molecules within 
the sphere recorded to calculate {p„}. Likewise, 
was constructed by monitoring the population within 
a coordination sphere of radius A around each water 
molecule. 

For coordination spheres with A > 3.5 A, po and xq are 
not reliably obtained directly from the simulation. Refer- 
ence values of these quantities were obtained by growing 
an empty cavity (for pq) or a coordination sphere (for xq) 
from A = 2.5 A to A = 5.5 A in steps of 0.25 A. For a cav- 
ity of radius A, the probability of the n = state within 
a shell of thickness 0.25 A was used to obtain the free 
energy to increase the radius from A to A + 0.25 A. This 
stratified sampling approach was used to construct the 
reference pq{X) and a;o(A) curves. (Ensemble reweighing 
can be used to obtain the probability in the shell but the 
results were not different.) 
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D. Reweighed ensemble simulations 

Each reweighed simulation was performed for 12x10^ 
sweeps of which the first 6x10^ were used for equilibra- 
tion. At the start of the equilibration phase all weights 
{rjn} were set to zero. After every j^^ set of 2000 sweeps, 
the occupancy statistics {pn} are obtained and a new set 
of weights for the (j + 1)*'^ sweep determined according 
to 

Vi"-' =vi- 0.1fcBTp„. (7) 

This cycle is continued till the end of the equilibration 
phase. (The factor of 0.1 serves to dampen the oscilla- 
tions in r]n.) In the subsequent 6x10^ sweeps of produc- 
tion, the set of weights obtained at the end of the equi- 
libration phase were used unchanged. The probabilities 
{Pn} are re-constructed from the reweighed simulations 
using Eq. [6j A similar strategy is followed for {a;„}. 

For large coordination spheres ( > 4.0 A ), because we 
have a large number of coordinating states, the reweigh- 
ing technique is adapted for overlapping subsets of {pn} 
and The entire distribution is then reconstructed 

by stitching together these subsets. (In Appendix I, 
we collect an explicit demonstration of the coordination 
states obtained using the reweighing technique together 
with a discussion of the thermodynamic meaning of the 
weights.) 



IV. RESULTS AND DISCUSSION 

The oxygen-oxygen radial distribution function reveals 
that four water molecules fill the first hydration shell 
(Fig. [2| . A closer examination of these four waters shows 
that there is asymmetry in their distribution about the 
distinguished water molecule. For example, the nearest 
water molecule is more localized than the fourth near- 
est water molecule. These differences also hint at differ- 
ent energetic interactions with the distinguished water 
molecule, an aspect that emerges below. 

Figure [3] shows the calculated hydration free energy 
of water for coordination radii between A = 2.5 A to 
A = 5.5 A, near the outer-edge of the second hydration 
shell. The agreement of the calculated over this broad 
a range of coordination radii serves as a rigorous check of 
the simulation. It also serves to ensure the correctness of 
the individual pieces — /ceTlnxo, kBTlnpo, and (0) 
— contributing to the excess chemical potential (Eq. |2]) . 
We next consider these individual pieces. 
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FIG. 2. Oxygen-oxygen radial distribution function (dashed 
line). The distribution of neighboring oxygen atoms is based 
on distance-order (filled circles) and their sum (open circles) 
is given as well. 
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FIG. 3. The hydration free energy of water using Eq. [2] for 
different coordination radii. The estimates are within ±0.35 
kcal/mol ~ ±0.6fcBr (shown by blue dotted lines) from the 
average estimate (black dashed line). kBTlnxo and kBTlnpo 
were obtained by calculating the work done in progressively 
enlarging the coordination sphere in steps of 0.25 A. The outer 
term is obtained using a 2-point Gauss-Legendre quadrature 
for electrostatic^^ and histogram overlaps* for van der Waals 
contributions. The solid black line is the estimate for A = A; 
a 3-point quadrature for electrostatics did not change this 
estimate by more than 0.2 kcal/mol. 



A. Outer-shell contribution /J^f (0) 

The outer-shell contribution (Fig. |4]), /i^(0), for var- 
ious coordination radii shows that for A > 3.5 A, the 
Gaussian model is in good agreement with the rigor- 
ous estimate based on Gauss-Legendre quadratures'^^ 



for electrostatics and histogram overlap!^ for van der 
Waals (together labeled as "Thermodynamic integra- 
tion" in Fig. |4]). When the distribution of binding en- 
ergy (e) of the solute water with the outer-shell material 
is GaussiarP^Jti^, the binding energy of the solute with 
the bulk is a sum of a large number of small, uncorre- 
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FIG. 4. Tlie outer-shell contribution to the hydration free 
energy (Eq.[2| for various coordination radii. Gaussian model: 
Result based on describing the distribution of solute-solvent 
binding energies as a Gaussian. Thermodynamic integration: 
Result based on Gauss-Legendre quadrature for electrostatics 
and histogram overlap for van der Waals interactions. 



lated contributions. Contrarily, the failure of the Gaus- 
sian model for A < 3.5 A indicates that water molecules 
at this range continue to experience specific interactions 
with the distinguished solute. Thus A w 3.5 A is the ra- 
dius to account for all the local, specific hydration effects. 
Interestingly, 3.5 A is also nearly the same cutoff distance 
used in analyzing hydrogen-bonding in liquid water (for 
example, Refs. [3 and EH). This agreement serves to 
emphasize the strength of the frameworlP^^tSl jji sorting 
specific, chemical and non-specific, bulk interactions in a 
physically transparent fashion. 



B. Packing and chemical contributions 

Figure [5] shows the packing and chemical contributions 
for various inner-shell radii. As compared to the revised 
scaled particle theory result!^, the packing estimates ob- 
tained here begin to increase faster with A for A > 4.5 A. 
This is expected, since the present simulations are at con- 
stant volume and system size effects will influence the 
opening of progressively large cavities. In particular, in 
a constant volume simulation, the spontaneous opening 
of a large cavity is expected to be inhibited relative to 
what would be observed in a constant pressure simula- 
tion. However, the same effects will enhance the chemical 
term as well, and thus a cancellation of errors can be ex- 
pected (Fig. [3]). 

Combining the packing and chemistry estimates, leads 
to two interesting observations. First, at A ~ 3.1 A, local 
chemical and packing contributions balance each other 
and the net hydration free energy is given solely by the 
outer-shell contributions. A similar behavior has also 
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FIG. 5. The packing (fcBTlnpo) and chemical contributions 
(kBTlnxo) for different coordination radii (A). The packing 
contribution predicted by the revised scaled particle theory^ 
is shown for comparison. po(A) and a;o(A) are obtained by a 
stratified sampling approach (see text). At 3.1 A the packing 
and chemical effects balance and the outer term, ^i™(0) ac- 
counts for the entire free energy of hydration. The deviation 
of packing estimates from the revised scaled particle theory 
result for A > 4.5 A indicates onset of system size effects in a 
constant volume simulation. 



been found for water simulated with the TIP3P water 
modefl^'^ and some ab initio models^^ 



Second, by 

about 5.5 A, the outer-shell contribution is small (Fig.|4|, 
and the difference in packing and chemical contributions 
itself accounts for nearly all of the hydration free energy 
of water. 



C. Molecular aufbau and water coordination 
structure 

Figure [6] shows the building-up of the local, chemi- 
cal contribution, fesTlnxo, with successive addition of 
a solvent ligand. At A = 3.5 A — the coordination ra- 
dius that fully accounts for the local, specific contribu- 
tion and beyond which outer-shell interactions are non- 
specific (Fig.|4| — approximately all of the local chemical 
contribution is accounted for by the n = 4 coordination 
state (Fig. |6j Left panel). Thus four water molecules 
around the distinguished water molecule are necessary 
to fully capture the local chemistry. However, the four 
solvent ligands are not equivalent. The first added lig- 
and, corresponding to the n — 1 coordination state, con- 
tributes about twice (thrice) as much to fee Tin xq as the 
second (fourth) ligand. The importance of the n = 1 co- 
ordination state is clearly seen for A up to 4.5 A (Fig. |6j 
Right panel). 

For coordination radii beyond A — 4.5 A, the distinc- 
tion between the solvent ligands tends to decrease, and 
by A = 5.5 A, both the n = 1 and the n = 4 states con- 
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FIG. 6. Left panel: The local chemical contribution with successive addition of a water molecule to the coordination sphere. 
ksTlnxo^i accounts for the chemical contribution upon adding i water molecules to the coordination sphere (Eq. [3|. The dotted 
line is the stratified sampling estimate from Fig. [5] The excellent agreement between these estimates serves as a rigorous check 
of the simulation data. The symbol size is about a fcaT and is a conservative estimate of the statistical uncertainty in the data. 
Right panel: The change in the ksT In xo upon changing the coordination state from i to i + 1. For clarity, only the — > 1 
and 3 — >■ 4 transitions in the coordination state are shown. 



tribute roughly equally to fceTlnxo (Fig.|6j Right Panel). 
(In contrast, such a turn-around is not expected in the 
hydration of ions-*^.) We next explore the reason for this 
behavior and discover the importance of ligand-bulk in- 
teractions in influencing the ligand's interaction with the 
distinguished solute. 

The maximum term approximation of Eq. [3] 

k^T \nxo i ~ — fcerin — + Wi 
Po 

shows that the intrinsic occupancy of the solute- free coor- 
dination sphere, the Pi/po term, and solute-ligand inter- 
actions, the Wi term, together determine the local chem- 
ical contribution. For coordination states n = 1, ... ,4, 
Pi/po > 1 a-nd hence the first term on the right will al- 
ways contribute favorably to the chemical contribution. 
Further, —kQThi{pi/pQ) is comparable to, and in some 
instances more favorable than, Wi itself (Fig. [t]) : thus in- 
trinsic occupancy variation, reflecting the bulk behavior 
of the fluid, is important in the hydration structure of a 
distinguished water molecule. 

Figure [t] (Left Panel) shows that Wi decreases for A < 
4.5 A. Observe that the free energy of forming a 1-water 
cluster is comparable to the free energy for forming a 4- 
water cluster. This observation reiterates the importance 
of the n = 1 coordination state in the thermodynamics 
of hydration (Fig. [6]). 

But for A > 4.5 A, W^ increases indicating that it is 
progressively less favorable to form a cluster around the 
distinguished solute. To understand why Wi increases 
for A > 4.5 A, let us consider —kBT\ii{pi/po) for var- 
ious A (Fig. [6j Right Panel). The curve for i = 1 is 
relatively flat, and for A > 4.5 A it is at a value of ap- 
proximately — 1.6 kcal/mol. This is a rather remarkable 



trend. The free energy is negative (favorable) indicat- 
ing that the solvent ligand within the empty coordina- 
tion sphere can satisfy the bonding of solvent molecules 
adjoining the coordination sphere. Indeed the near con- 
stancy of the value suggests that the solvent ligand in 
the coordination sphere is near the interface with the 
outer-shell solvent. This fact is easily confirmed by cal- 
culating the pair-correlation of the solvent ligand and 
remaining water molecules (data not shown). The mag- 
nitude of the free energy is consistent with that expected 
for hydrogen bonds^ and is nearly 20% of the hydra- 
tion free energy of SPC/E water. This observation com- 
plements one made recently by Chempath and PratlP^ 
that the bulk of the hydration free energy of water is ac- 
counted for before the water molecule fully enters the liq- 
uid from the vapor. Even more remarkably, for A > 4.5 A, 
— /cb 2^111(^1 j/po) « —1.6 X i; that is Pi oc p\. Thus for 
A > 4.5 A, the solvent ligands entering the solute-free co- 
ordination sphere are largely independent of each other 
and stay near the interface. 



The above discussion then provides an answer to why 
beyond A > 4.5 A the aufbau approach finds that the 
distinction between the contributions of the coordinat- 
ing solvent ligands decreases (Fig. [6j Right Panel). As 
the coordination sphere is enlarged, the solvent ligands 
prefer to partition to the interface. That is the solute 
water-solvent ligand interaction is not strong enough to 
compete with the bonding opportunities for the ligand 
with the outer-shell material. 
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FIG. 7. Left panel: Free energy to form an i-water cluster around a distinguished water molecule in a coordination sphere 
of radius A. Right panel: Intrinsic occupancy variation in the coordination sphere for various coordination states and 
coordination radii. 



V. CONCLUDING DISCUSSION 

The number of solvent ligands within a coordination 
sphere around the solute, the coordination state of the 
solute, depends on two factors. One is the ligand oc- 
cupancy of the coordination sphere in the absence of the 
solute and the second is the free energy of forming solute- 
ligand clusters within the coordination sphere. For the 
hydration of a water molecule solute, these two factors 
are of comparable magnitude. 

For a chemically meaningful coordination radius of 
3.5 A, the free energy of forming a 1-water cluster around 
a solute water is nearly the same as that for forming a 4- 
water cluster, emphasizing the importance of the n — I 
coordination state in the hydration thermodynamics of 
the solute water. Yet the average coordination state of 
the solute water is n = 4. We can understand this behav- 
ior by noting that we incur an energetic penalty to create 
an empty coordination sphere in bulk water. Solvent lig- 
ands that enter this empty volume partition towards the 
surface of the coordination sphere, thereby satisfying the 
bonding requirements of the water molecules outside the 
coordination sphere and also shrinking the coordination 
volume itself. This effect itself accounts for nearly 2-3 
water molecules of the 4 water molecules comprising the 
average coordination structure of the solute water. It is 
precisely because the bonding between the solute and the 
solvent ligand is strong enough to only organize a 1-water 
cluster that the solvent ligands also bond with the outer- 
material, where more bonding opportunities are avail- 
able. 

These results can be better appreciated by consider- 
ing the hydration of strongly ionic solutes. There the 
local ion-solvent water interaction overwhelms the inter- 
action of the solvent ligand with the outer-material. Thus 
the occupancy distribution of the empty coordination 



sphere plays but a minor role in the observed coordina- 
tion structur^^. For example, for a coordination radius 
of 3.5 A, the free energy for forming a one-water cluster 
around Na+ is nearly 10 times the contribution from the 
intrinsic occupancy variation, — fceTlnpi/po- For wa- 
ter, in contrast, Wi is nearly the same as — fceT'lnpi/po- 
Thus the hydration of a solute water molecule results 
from an interplay of two small, energetically equivalent 
contributions: one being the direct interaction of the dis- 
tinguished water molecule with the proximal solvent lig- 
and, and the second being the interaction of the proximal 
solvent ligand with the bulk. 

Our results suggest that we cannot claim that the aver- 
age coordination structure is not four because the free en- 
ergy of forming a 1-water cluster is comparable to that for 
forming a 4-water cluster around a solute. That is, given 
that the bonding of one water to the solute is as strong as 
the bonding of four does not invalidate the prominence 
of four coordinate states in the bulk fluid. This is espe- 
cially so when solute-ligand interaction is comparable to 
ligand-bulk interaction as happens for water modeled by 
SPC/E. 

The spectroscopy experiments (on real water j^EI ap- 
pear to probe the local bonding of water. Our results 
on SPC/E could possibly prove helpful in better appre- 
ciating those results. The development presented here 
together with advances in modeling the spectrgP^for sim- 
ulation models of water could also prove fruitful in better 
appreciating liquid water. 



ACKNOWLEDGMENTS 

We thank Purushottam Dixit for helpful discussions. 
D.A. thanks the donors of the American Chemical Soci- 
ety Petroleum Research Fund for the financial support. 



8 



Appendix A 



Below we present an explicit example of reconstructing 
the coordination states based on the ensemble reweighing 
technique (Section II A). We also establish the thermo- 
dynamic meaning of the weights (Eq. [5| . 

Figure [S] shows the raw and reconstructed {a;„} dis- 
tribution for an inner-shell of radius 4.0 A. Similar re- 
sults are obtained for {pn} and for different coordination 
radii. As Fig. [8] shows, the reweighing approach helps 
uncover the low probability states. Comparing the raw 
{xn} and reweighed {a;„} distributions, shows that the 
weights {?7j,w} increase the instances of rarely occurring 
states at the expense of the more probable ones. Note 
that in the reweighed simulation, the increased fluctua- 
tion of population in the coordination sphere can degrade 
the accuracy in estimating the probability of states, but 
together with the computed weights, a robust estimate 
of the probabilities {a;„} can be formed. 
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FIG. 8. The distribution of coordination states, {xn}, around 
a distinguished water molecule. The coordination radius A is 
4 A. From the raw populations obtainable from the simula- 
tion record, the n — state was not observed. The reweighed 
probabilities, {x„} (blue squares) are close to the desired 
equiprobable distribution. Reweighing (red circles) uncovers 
xo- 



Besides serving as aides in calculating the probabilities, 
the weights have a physical meaning. If the weights are 
chosen such that the {pj} and {xj} become uniform, then 



Eq. [6] implies that 



Pj 



Pr& 



Xje^"^''^ 
x^pP-n^-^ 



= 1 



Xj/pj 



o0(')j-'?j,w) 



Xr/Pr e^('''— 

But from Eq. [l| we have 

Xj/p 
Xr/P' 
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FIG. 9. Relationship between the weights {»7j} and {t/j.w} 
and the excess chemical potential in a defined coordination 
state, /i^(n), for a coordination sphere of radius A = 4 A. The 
most probable state n = 8 (Fig. [8| is taken as the reference 
(r). The difference of weights in reweighed simulations (red 
squares) agree reasonably with the difference in hydration free 
energy of the solute in the corresponding coordinate states 
(black open squares). The agreement is not perfect since the 
requirement Eq. |Al| is not exactly satisfied and the explicitly 
calculated values /i™(n) also have uncertainties of about a 



Comparing Eq. [A2] and Eq. |A3j 

/^w (i) - (^) = ('7j,w - %) - (?yr,w - Vr) , (A4) 

where r is the selected reference state. Thus the differ- 
ence of the weights for the {p„} and {xn} distributions 
for any coordinate state j equals the excess chemical po- 
tential of the solute in that defined coordination state 
to within a constant dependent on the chosen reference 
state (Figure l9|. 
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